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ABSTRACT 


This j^per outlines the equations which govern the 
non isothermal flow of reactive fluids, and describes the means by 
which finite element analyfiis can be used to solve these 
equations for th% sort of arbitrary boundary conditions encount- 
ered in industrial practice. The performance of the computer 
code is illustrated by several trial problems, selected more for 
their value in providing insight to polymer processing flowr than 
as practical production problems. Although a good deal remains 
to be learned as to the performance and proper use of this numer- 
ical technique, it is undeniably useful in providing better 
understanding of today's complicated polymer processing problems. 


INTRODUCTION 
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Finite element analytis- of fers great promiie for reducing the 
empiricism often found in polymer processing design, since it is 
well suited for modeling the complicated boundary conditions and 
material properties encountered in industrictl practice. Although 
the method is now well accepted in structural stress analysis, 
its use in fluid transport situations is less widespread. Ne 
have sought to exploit some of the many advantages of the finite 
element method in polymer fluid processing analysis, and this 
paper describes some of our work in chemorheology . 

Reactive flows have several advantages in polymer processing 
in comparison with more traditional melt'flow techniques. Per- 
haps the most significant of these is the energy savings which 
result from the elimination of the several melting stages neces- 
sary in such technologies as extrusion or injection molding. 
However, these advantages are offset to some degree by the com- 
plexity of the process, which renders the empirical approach to 
process development difficult in the extreme. The flow veloci- 
ties in such processes are governed by the fluid viscosity, which 
is a strong function of temperature and molecular weight. The 
temperature in turn is affected by the viscous dissipation and 
the heat released or consumed by the reaction, and the reaction 
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rate is also a strong function of temperature. All of these var- 
iables interact and change in such a way as to make an intuitive 
grasp of the process almost impossible, and there is obviously an 
advantage to being able to provide some sort of mathematical or 
numerical simulation of the process. 


The finite element scheme to be described below is very useful 
for obtaining numerical solutions for reactive flows with arbi- 
trary boundaries, and such a technique is well suited for 
detailed analysis of real industrial processes. However, we 
argue that the greatest value of the method may not be in these 
detailed calculations, but in the degree to which the method can 
enhance the process designer's intuition. For this latter pur- 
pose, it is often sufficient to run only very small and 
inexpensive trial problems; these can elucidate the manner in 
which the various problem parameters interact, so that the 
designer develops a much improved "feel" for the problem. It is 
likely that most real advances in processing technology will con- 
tinue to be made more by Edisonian innovation than by detailed 
mathematical calculations. But since today’s processes have 
become so complicated, « such a technique as the finite element 
method can be a powerful adjunct to intelligent intuition. 


ORIGINAL PAGE IS 
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The equations which govern the non iso thermal flow of a reac* 
tive fluid are derived in several texts on transport phenomena 
and polymer processing (e.g. refs. 1,2). Regarding velocity, 
temperature, and concentration of unreacted species as the f undr.- 
mental variables, the governing equations can be written as: 

^[au/at + u«^u] ■ -^p + »^*u 

^c[aT/at + u«^T] ■ Q + lc'^*T 

[aC/at ♦ u*^C] - R + D^*C 

(A list of nomenclature is attached.) The similarity of these 
equations is clear. In all cases, the time rate of change of the 
transported variable (velocity, temperature, or concentration) is 
balanced by the convective or flow transport terms (e.g. u«^C), 
the diffusive transport (e.g. 0^*C), and a generation term (e.g. 
R). 


The analyst seeks expressions for the space- and time-depen- 
dent velocities, temperatures and concentrations which satisfy 
these equations and also the problem's boundary conditions. It 
is generally the boundary conditions which make real problems 
intractable: even if one were able to describe the boundaries 
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mathcfflatically , th« result ing sxprtssions would not liksly bt 
snsnablc to closad~form solution. In addition, many of the "con- 
stants" in the above equations are often nonlinear functions of 
the problem variables. In reactive polymer processing, one might 
encounter such expressions as the following: 


n - to»"“’«*PCBi/R 9 T]exp{/p)(mw» <) 

Q - (n/2)i%t%) ♦ R(AH) 

R • -k,exp(Ba/RgT)C» 


Zt is clear that all of these expressions, taken together, 
constitute a mathematical situation which must be approached with 
caution. Even though it is not overly difficult to incorporate 
them into a numerical scheme, which we have done, it is important 
to proceed slowly enough to develop the proper experience in the 
computer code's behavior before tackling full-blown problems. 

COMPOTER MODEL 

We have sought to develop a finite element code which is able 
to predict polymer fluid velocities, stresses, temperatures, and 
degrees of chemical conversion in a variety of flow geometries 
and for a variety of fluid material properties. Space limita- 
tions prohibit our listing here the full derivation of the finite 
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tl«fflent «qu«tioni from tht abovt diffarantial tquationi, but 
thara ax 1st lavaral wall davalopad naans by which this nay ba 
dona. Tha raadar is diractad to standard taxts (3,4) for a nora 
conplata dascription, and wa will just stata. hare that wa anploy 
tha Galarkin nathod of weighted residuals, together with 
isoparanatric napping and interpolation, to replace tha differan* 
tials by integrals which can ba evaluated nunarically over nnall 
subregions ("alanants”) . Tha results of these numerical inte- 
grations are than assanblad into a sat of simultaneous algebraic 
aquations which can ba solved nunarically. 

Tha salient features of our coda can ba listed briefly as fol- 
lows: (1) Velocity, tamparatura, and chemical conversion are 
taken as nodal unknowns, so that coupled incompressible viscous 
flow and diffusive-convective heat and mass transport may ba mod- 
eled. (2) Incompressibility is enforced by a "penalty* formu- 
lation employing selective reduced integration. This approach 
requires the use of double precision arithmetic with a 
concommitant reduction in the amount of available core, but it 
has several programming advantages which usually outweigh this 
drawback. (3) The code is developed primarily for plane and 
axisymmetric flows. We have coded a three-dimensional 
capability, but generally feel that the expense of running 
three-dimensional problems is not justifieo for most of our mod- 
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•ling retearch. (4) The code includes several models for the 
effect of shear rate, temperature, pressure, and chemical convert 
sion on the fluid viscosity. These nonlinear models have not yet 
been researched extensively, hovever, and the exploration of the 
interative schemes needed for their proper use constitutes a 
major goal for future research. We have also coded a capability 
for viscoelastic fluid effects (ref. 5), but currently feel that 
this difficult aspect of polymer flow is being researched satiS' 
factorily by other workers. (5) Convective transport of heat or 
chemical species can be handled either by a conventional Galerkin 
treatment or by the convenient but still controversial "optimal 
upwinding” approach. (6) The code can treat transient problems 
by means a two-point "theta-method” time-stepping algorithm. The 
dynamic -algorithm is also useful in nonlinear problems, in which 
the final fluid state may be approached dynamically from an esti- 
mated initial state. (7) The code is capable of a variety of 
iterative treatments of nonlinear problems, including 
Newton-Raphson iteration and incremental load methods. 


Some additional discussion is warranted concerning the treat- 
ment of convective effects beyond what has been mentioned in item 
(5) above. Momentum convection (^u*^u) is generally negligible 
in comparison with the viscous terms due to the high viscosities 
of polymer fluids, but the convective terms tend to dominate the 
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energy and mass transport equations due to the generally low 
thermal and mass diffusivities* The prograouning of the 
convective terms presents, no special problems in the Galerkin 
approach beyond the need to store and solve unsymmetric matrices, 
but it is well known that the presence of strong convective terms 
tends to create oscillations in the final solutions which can be 
large enough to destroy their value. This instability is related 
to the tendency of convection to produce large downstream gradi- 
ents which the finite element grid cannot resolve. A largely ad 
hoc procedure known as "upwinding" has been used in both finite 
element and finite difference work which seems to alleviate this 
problem by providing a greater numerical weight to the upstream 
portion of the element. Hughes (6) has published a very conven- 
ient means of upwinding, in which the sampling points in the 
numerical integration scheme are simply moved upstream an appro- 
priate distance. We have made extensive use of the Hughes 
upwinding technique, but the reader is cautioned that this method 
is regarded as controversial by many workers. A provocative 
paper by Gresho (7) details many of the possible pitfalls in 
upwinding, and states a strong preference for grid refinement as 
the appropriate cure for convection-induced instabilities. 

SELECTED NUMERICAL RESULTS 


9 


PAGE a 

Of POOR QUALITY 


Nonreactive entry flow. Figure 1 shows the streamlines for a 
4tl entry flow. Here a grid of 100 four*‘node linear elements was 
used to model the upper symmetric half of a plane capillary, and 
a fully-developed Poiseuille velocity was imposed on the reser- 
voir entry as a boundary condition. The streamlines are identi- 
cal with published experimental and numerical results, although 
the grid used here was not intended to be fine enough to capture 
the weak recirculation which develops in the stagnant corner of 
the reservoir. 


The temperature contours for convectionless flow are shown in 
figure 2, which shows a hot region at the entrance of the 
capillary due to the combination of high viscous energy dissi- 
pation there and its distance from cool boundaries to which heat 
may be conducted. These isotherms are normalized on the maximum 
centerline temperature expected for Poiseuille flow in the 
capillary. 

The importance of thermal convection relative to conduction is 
given approximately by the Peclet number Pe ■ UL^c/k, where U and 

ft 

L are a characteristic velocity and length. Figure 3 plots the 
variation of temperature along the centerline for various values 
of the Peclet number, and it can be seen that the effect of 
increased thermal convection is to sweep the cooler upstream flow 
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particlts into th« capllltry, with • rotulting lowering of the 
temperaturei overall and a ahift downatrean of the hot apot near 
the throat. The relatively coarae grid uaed in thia problem 
produced unatable Galerkin reaulta for Peclet numbera higher than 
approximately ten, and the higher degreea of thermal convection 
were computed uaing the upwinding formulation. Further teata 
with refined grida ahould be completed to aaaeaa the accuracy of 
the upwinded aolutiona, although the plota in figure 3 appear 
reaaonable. 

One-dimenaional reactive flow. Aa a preliminary trial problem 
in our computationa of reactive flow, we have atudied a aimple 
aituation in which a fluid obeying firat-order chemical kinetica 
movea at conatant velocity and temperature in the poaitive x 
direction. Here only the maaa-tranaport equation- ia operative, 
and it takea the aimple formt 

u(dC/dx) - -KC ♦ D(d‘C/dxM 

Thia equation ia aolved eaaily, and for nonaeto valuea of the 
diffuaion coefficient D two boundary valuea for C muat be aped- 
fied. One of theae ia the initial concentration at the inlet, 
and the other requirea a conaideration of the outlet conditiona. 
Here aeveral poaaibilitiea exiat, and we have atudied the caae in 
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which th« conetntration of tha outlet reiervoir it allowed to 
rite to meet that tupplied by the flow; thit it equivalent to 
tpecifying a zero concentration gradient at the outlet. For the 
cate of negligible diffution the tecond**order term van- 
ithet from the above equation and the downttream boundary condi- 
tion cannot be tpecified. The aolution it then a aimple 
exponential, in which the reactive tpeciea vanith according to 
firtt-order kinetica at they are carried downttream at conttant 
velocity. 


Figure*4 thowt the computed and exact predictiont for reactive 
group concentration at a function of dittance along the channel. 
The Galerkin valuea are nearly exact, but it it clear that 
upwinding leadt to er«roneout retultt in the tmall-dif fution cate. 
The upwinding hat introduced an artificially high diffutivity and 
a zero concentration gradient at the outlet, and tuch artifacta 
mutt be conaidered at potaibilitiea when upwinding it uted. 

Two-dimentional nonitothermal reactive flow. Figure 5 thowt 
the contourt of conttant converaion for a tvo-dimentional analog 
of the flow diacutted in the previout tection. Here again, a 
timple uncoupled problem it treated in which the material parame- 
tera are taken to be independent of the aolution variablet, and 
in which the velocity conditiont are pretcribed. The concen- 
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tration ii takan to havt a fixtd valut at tht inltt and a zaro 
gradiant at tha outlat, aa bafora. Tha two-dimanaionality of tha 
problan ii eontainad In two faaturaa: tha valocity la takan to 
ba parabolic r ranging from a maximum at tha cantarlina to zarc at 
tha walla (a Poiaauilla flow); and now diffiaiva haat and maaa 
tranaport can occur in both tha x and y diraction. For tha low 
dif fuaivitiaa ahown, maaa diffusion in tha x diraction is nagli- 
gibla, as was daitfonatrar.id in tha previous aaction. Howavar, tha 
concantration gradiants in tha y diraction ara substantial, so 
that diffuaiva transport in that diraction is appraciabla avan at 
0 • 0.01. At D ■ 0.001, avan tha y-diffusion ia nagligibla, so 
tha concantration contours simply raprasant a fluid which movas 
in tha x-diraction whila raacting by first-ordar kinatics. Tha 
concantrations along tha cantarlina art than idantical with tha 
D ■ 0.01 curva of figura 4. 

Pigura 6 shows tha contours of constant tffflpara^ure which 
rasult from this flow (with D ■ .001), whara tha tamparatura 
boundary conditions wars sat to zaro at tha antry and along tha 
top and bottom nurfatas. Tha tamparatura gradiant at tha outlat 
was allowad to become zaro, similar to tha conrentration 
gradiant. Tha rasults obtainad for tha tamparatura fiald ara of 
coursa dapandart on tha valuas chosan for fluid propart ias. To 
avoid using spaca hara for a datailad discussion of tna diman> 
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sional analysis used for selecting these parameters, we will 
state simply that in figure 6 the viscous dissipation and 
reaction heat make approximately equal contributions to the 
internal heating (Brinkman and Damkoehler numbers both equal to 
three). 
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CONCLUSIONS 


The numerical model described above has a significant present 
ability to simulate a wide range of problems in polymer process- 
ing. At the same time, it is small enough to permit easy imple- 
mentation in even rather small processing facilities, and for 
quick familiarization by process designers. We feel such a capa- 
bility would lead to a significant advance in industry capability 
for process development and optimization. 
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NOMENCLATURE 


c specific heat 

C concentration of reactive species 
D diffusion coefficient 
El activation energy for viscosity 
Ea activation energy for reaction 
Ic coefficient of thermal conduction 
km reaction rate preezponential factor 
R overall reaction rate 
m reaction order 
mw molecular veight 
n power *’law exponent for viscosity 

p pressure 

Q Internal heat generation 
R internal species generation 
Rg gas constant 
T temperature 

u velocity or velocity vector 
A factor for pressure dependency 
% shear rate 
f viscosity 
Vo viscosity coefficient 
p density 
AH heat of reaction 
^ Laplacian operator 
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Figure 6 - Isothermal contours in two-dimensional reactive flow, 
Galerkin calculation. Heat generation by viscous dissipation and 
reaction heating. 



